Journal of Chemical Theory and Computation

Accurate ion force fields are essential for molecular dynamics simulations of biomolecular systems, particularly in combination with modern water models such as OPC. While OPC water improves the description of bulk water and biomolecules, the transferability of existing ion force fields to this model remains an open question. Here, we systematically assess the transferability of monovalent and divalent ion force field parameters (Li+, Na+, K+, Cs+, Mg2+,Ca2+, Sr2+, Ba2+, Cl- and Br-) to OPC water by comparing single-ion and ion-pairing properties with experimental data. Our analysis reveals that no single literature parameter set provides accurate results for all ions when directly transferred to OPC water. We hence introduce the MS/G-LB(OPC) force field, which combines Mamatkulov-Schwierz-Grotz cation parameters with Loche-Bonthuis anion parameters. MS/G-LB(OPC) reproduces hydration free energies, first-shell structural properties and activity derivatives at low salt concentrations. Our results demonstrate that transferring ion parameters to OPC can lead to significant and ion-specific deviations from experimental data, making careful validation essential. At the same time, the systematic transfer and combination of ion parameters from existing force fields can provide a practical and computationally efficient alternative to full reparameterization. MS/G-LB(OPC) is available at https://git.rz.uni-augsburg.de/cbio-gitpub/opc-ion-force-fields.

Understanding RNA structural dynamics is essential for elucidating its biological functions, and molecular dynamics (MD) simulations provide an important atomistic complement to experimental approaches. However, the predictive power of MD is fundamentally limited by the accuracy of the underlying empirical Force Fields (FFs), particularly in capturing the delicate balance of non-bonded interactions. Here, we present a systematic reparameterization strategy that replaces the external gHBfix19 hydrogen-bond (H-bond) correction potential with an equivalent set of NBfix Lennard-Jones modifications within a state-of-the-art RNA FF. Using a quantitatively converged temperature replica-exchange MD ensemble of the GAGA tetraloop, we employed a reweighting-based optimization protocol to derive NBfix parameters that reproduce the thermodynamic effects of the original gHBfix19 terms. Sequential optimization of individual gHBfix19 components proved essential to ensure stable and transferable parameter refinement. The resulting fully reformulated NBfix-based variant, termed OL3CP-NBfix19, was validated on a representative set of RNA motifs, including tetranucleotides, A-form duplexes, and tetraloops. Across all tested systems, its performance is comparable to that of the reference gHBfix19 FF. By embedding the H-bond corrections directly into the standard non-bonded framework, the NBfix formulation eliminates external biasing potentials, simplifies practical deployment, and reduces computational overhead. Beyond this specific reparameterization, our results demonstrate a practical workflow for translating targeted H-bond corrections into native FF terms for efficient biomolecular simulations.

Molecular dynamics (MD) simulations are a powerful tool for investigating biomolecular dynamics underlying biological functions. However, the accessible spatiotemporal scales of conventional all-atom simulations remain limited by high computational costs. Coarse-graining reduces these costs by decreasing the number of interaction sites and enabling longer timesteps. In extreme cases, proteins are represented as single spherical particles; while such approximations facilitate cellular-scale simulations, they often sacrifice essential structural information, such as molecular shape and interaction anisotropy. Here, we present CGRig, a rigid-body protein model with residue-level interaction sites designed for long-time, large-scale simulations. In CGRig, each protein is treated as a single rigid-body embedding residue-level interaction sites. Its translational and rotational motions are described by the overdamped Langevin equation incorporating a shape-dependent friction matrix. Intermolecular interactions are calculated using G[o]-like native contact potentials, Debye-Huckel electrostatics, and volume exclusion. We validated that CGRig accurately reproduces the translational and rotational diffusion coefficients expected from the friction matrix for an isolated protein. For dimeric systems, the model successfully maintained native complex structures. Furthermore, two initially separated proteins converged into the correct complex with an association rate consistent with all-atom simulations. Notably, CGRig achieved a simulation performance exceeding 17 s/day for a 1,024-molecule system. These results demonstrate that CGRig provides an efficient framework for simulating protein assembly while retaining residue-level interaction specificity, making it a valuable tool for investigating large-scale biomolecular self-assembly.

Protein kinases regulate signaling by recognizing short sequence motifs, and how these motifs bind influences both specificity and therapeutic strategies that target kinase pathways. Peptide-based inhibitors that engage substrate-recognition regions are attracting interest, but designing them requires an understanding of how a flexible peptide approaches and settles into the bound pose. Traditional studies have focused on the bound pose and affinities, whereas the steps that link the initial encounter with the bound pose have been explored less thoroughly because the relevant intermediates are too short-lived to capture experimentally and evolve on timescales that standard molecular dynamics cannot readily access. Here, we focused on Abl kinase and Abltide, the experimentally identified optimal substrate peptide for Abl kinase, and examined the sequence of events linking initial encounter to the bound pose using two-dimensional replica exchange (gREST/REUS), which selectively enhances flexibility in the peptide and its binding interface while also sampling progression along a distance coordinate. The resulting simulations yielded a detailed binding landscape, revealing five distinct encounter regions outside the substrate-binding site and six intermediate states that may connect the initial approach to the bound pose. Some encounter regions and intermediate states participate in the dominant binding pathways. During this process, EF/G/{beta}11 hydrophobic patch, together with G helix negative patch, plays a central role in guiding Abltide toward the substrate-binding site. These findings provide mechanistic insight into substrate recognition by protein kinases and offer a foundation for the rational design of peptide-based inhibitors.

Protein function often involves multiple conformational states. Several multiple sequence alignment-perturbing strategies, including stochastic subsampling, clustering, and column masking, have been shown to enhance AlphaFold2 (AF2) sampling of alternative protein states. Here, we evaluate these strategies on AlphaFold3 (AF3) and compare their performance with the BioEmu Boltzmann sampling model on 107 proteins with multiple experimentally solved conformational states. We find that unperturbed AF3 samples alternative states with significantly higher TM-scores compared to AF2 and comparable to BioEmu. In particular, all MSA perturbation methods improve AF3 sampling at a statistically significant level, improving the top 1% TM-score by at least 0.05 in approximately 20% of cases each, while rarely worsening the performance. Furthermore, we find that different choices of amino acid masks can improve column-masked AF3 sampling for specific targets. Our results highlight how MSA perturbations remain relevant in AF3, providing a useful tool for understanding dynamic biological processes.

The rational design of biomolecule immobilization strategies requires molecular-level understanding of how surface properties, tethering geometry, and structural dynamics jointly influence stability and function. Recently, coarse-grained molecular dynamics simulations based on the Martini force field have emerged as an efficient framework for studying enzyme-surface interactions. However, the reproducible construction of immobilized systems with controlled orientations remains technically challenging, usually involving multiple computational tools. Here we present MartiniSurf, an open-source command-line framework for the preparation of protein and DNA systems immobilized on solid supports within the Martini paradigm. MartiniSurf integrates automated structure retrieval and cleaning, coarse graining via tools from the Martini force field software ecosystem, customizable surface generation, and biomolecule orientation based on user-defined anchoring residues, producing complete GROMACS-ready simulation systems. The framework supports both implicit restraint-based anchoring and explicit linker-mediated immobilization, including surfaces functionalized with user-defined ligands or linker-like moieties, enabling representation of mono- and multivalent attachment geometries at different modeling resolutions. Structure-based G[o]Martini potentials can be incorporated for proteins, while DNA systems are modeled using Martini 2. Optional substrate insertion, pre-coarse-grained complex handling, and automated solvation and ionization further extend system flexibility. By integrating these components into a unified workflow, MartiniSurf enables systematic and high-throughput in silico exploration of surface-tethered biomolecules and provides a robust computational platform for rational immobilization studies. TOC Graphic O_FIG O_LINKSMALLFIG WIDTH=200 HEIGHT=146 SRC="FIGDIR/small/714767v1_ufig1.gif" ALT="Figure 1"> View larger version (45K): org.highwire.dtl.DTLVardef@bc1ac4org.highwire.dtl.DTLVardef@1813b43org.highwire.dtl.DTLVardef@159b19borg.highwire.dtl.DTLVardef@19b60d6_HPS_FORMAT_FIGEXP M_FIG C_FIG

There is a growing demand for molecular dynamics simulations to explore longer timescale behavior of giant protein-DNA complexes such as chromatin. To address this need, we extended OpenCafeMol, a GPU-accelerated residue-level coarse-grained molecular dynamics simulator originally developed for proteins and lipids, to support 3SPN.2 and 3SPN.2C DNA models. We also implemented a hydrogen-bond-type many-body potential to model DNA-protein interactions more accurately. To further improve computational efficiency, we introduced a localized scheme for calculating base-pairing and cross-stacking interactions. Benchmark tests show that OpenCafeMol on a single GPU achieves up to 200-fold speed-up for DNA-only systems and up to 100-fold speed-up for DNA-protein complexes compared to CPU-based simulations. To demonstrate the capability of our implementation for long-timescale biological processes, we simulated an archaeal SMC-ScpA complex undergoing DNA translocation via segment capture (a proposed mechanism for DNA loop extrusion) in the presence of a DNA-bound obstacle. We observed continuous captured-loop growth accompanied by obstacle bypass within the segment capture framework.

Recently, a number of tools have been released that generate ensembles of protein structures based on artificial intelligence (AI) approaches. Although ensembles generated by the tools differ significantly, we demonstrate a computational path to harmonizing the various outputs under a stationary condition using two complementary physics-based approaches. In the first stage, the AI ensemble is used to seed a weighted ensemble (WE) simulation, promoting relaxation toward the steady state. In the second stage, trajectory segments generated by WE are reweighted to steady state using the recently developed RiteWeight (RW) algorithm. We applied this approach to generate an atomically-detailed equilibrium ensemble of unliganded adenylate kinase conformations, starting from ensembles produced by three AI tools: AFSample2, ESMFlow-PDB (trained from PDB structures), and ESMFlow-MD (trained from molecular dynamics simulation data). Dramatic differences in the AI-generated ensembles are largely erased during the WE-RW process, yielding a consistent description of the equilibrium ensemble for a given force field.

Biomolecular condensates formed by complex coacervation of highly charged proteins provide a powerful framework to understand how microscopic interactions give rise to macroscopic material properties. Atomistic molecular dynamics simulations provide detailed insights but remain limited in accesing the spatio-temporal scales relevant for condensate behavior. Here, we use the residue-level coarse-grained Mpipi-Recharged model to investigate condensates formed by ProT and positively charged partners, including histone H1, protamine, poly-lysine, and poly-arginine. Material properties, in this context, provide a stringent experimental benchamark for coarse-grained models. Our model reproduces salt-dependent phase behavior, protein binding affinities, and sequence-specific stability trends in agreement with in vitro experiments, despite the fact that material properties were not included in the model parametrization. We then establish a direct link between protein dynamics and macroscopic material properties by quantifying monomeric diffusion, conformational reconfiguration, and translational mobility within the dense phase, and relating these to condensate viscosity. By comparing dynamics across dense and dilute phases, we uncover a pronounced length scale-dependent behavior. While residue-level binding and unbinding events remain equally fast in both phases, protein reconfiguration time and self-diffusion are significantly slowed down within the condensates. This decoupling reveals how fast intermolecular interactions coexist with slow mesoscale condensate dynamics depending on the molecular length scale. Together, our results establish a predictive framework that links encoded sequence intermolecular forces and multiscale dynamics to the emergent material properties of complex biomolecular condensates.

The apparent pKa of ionizable lipids in lipid nanoparticles (LNPs) is a key determinant of RNA encapsulation during formulation and endosomal release after cellular uptake. However, it is difficult to predict the effective pKa of a given ionizable lipid solely from its solution pKa, because it is sensitive to the membranes composition, as well as solution conditions such as the salt concentration. We developed a simple continuum electrostatics model, based on Gouy-Chapman theory, to predict the shift in effective pKa for ionizable lipids in lipid bilayers as a function of salt concentration and membrane composition. We derive equations for the surface potential and fraction of lipids charged, which are solved self-consistently as a function of solution pH to extract the titration curve and effective pKa. The model shows that the shift in effective pKa is largest when the concentration of titratable lipid is high, and the effect is diminished by increasing salt concentration. We provide a python implementation of the model and an interactive notebook that will allow users to further easily explore the predicted pKa shifts as a function of formulation variables.

AI-based molecular dynamics simulation brings ab initio calculations to biomolecules in an efficient way, in which the machine learning force field (MLFF) locates at the central position by accurately predicting the molecular energies and forces. Most existing MLFFs assume localized interatomic interactions, limiting their ability to accurately model non-local interactions, which are crucial in biomolecular dynamics. In this study, we introduce ViSNet-PIMA, which efficiently learns non-local interactions by physics-informed multipole aggregator (PIMA) and accurately encodes molecular geometric information. ViSNet-PIMA outperforms all state-of-the-art MLFFs for energy and force predictions of different kinds of biomolecules and various conformations on MD22 and AIMD-Chig datasets, while adapting the PIMA blocks into other MLFFs further achieves 55.1% performance gains, demonstrating the superiority of ViSNet-PIMA and the universality of the model design. Furthermore, we propose AI2BMD-PIMA to incorporate ViSNet-PIMA into AI2BMD simulation program by introducing "Transfer Learning-Pretraining-Finetuning" scheme and replacing molecular mechanics-based non-local calculations among protein fragments with ViSNet-PIMA, which reduces AI2BMDs energy and force calculation errors by more than 50% for different protein conformations and protein folding and unfolding processes. ViSNet-PIMA advances ab initio calculation for the entire biomolecules, amplifying the application values of AI-based molecular dynamics simulations and property calculations in biochemical research.

Phase separated condensates are recognized as a ubiquitous mechanism of spatial organization in cell biology. Biophysical modeling of condensates provides critical insights into the dynamics and functions of these subcellular structures that are difficult to extract via experiments. Here we present an efficient computational pipeline, CASPULE (Condensate Analysis of Sticker Spacer Polymers Using the LAMMPS Engine), to simulate and analyze the biological condensates made of sticker-spacer polymers. CASPULE implements a unique force field that combines traditional Langevin dynamics with a "detailed balance proof" protocol for single-valent bond formation between stickers. This framework allows us to study the non-trivial biophysics that emerge out of the single-valent sticker interactions coupled with the effect of separation in energetic contribution by stickers and spacers. We provide detailed documentation on how to setup the simulation environment, perform simulations and analyze the results. Through case studies, we highlight the utility and efficacy of our pipeline. Importantly, we provide statistical parameters to characterize the cluster size distribution often observed in biological systems. We envision this tool to be broadly useful in decoding the interplay of kinetics and thermodynamics underlying the formation and function of biological condensates.

Intercalation of small molecules between DNA base pairs affects DNA conformation, disrupting essential cellular processes including replication, transcription, and repair. We investigated conformational changes in 18-mer DNA upon intercalation of doxorubicin, SYBR Gold and YOYO-1 using extensive MD simulations. Two main patterns for the intercalation were identified: RISE-type intercalation occurs between adjacent base pairs and extends the DNA helix with decreased twist angles, while OPEN-type intercalation proceeds through base-pair opening without significant DNA extension. Kinetic analysis revealed that association rates for intercalation followed the order: first YO-moiety (mono-intercalation) > SYBR Gold > doxorubicin > YOYO-1 (bis-intercalation). Free energy landscape showed that forces at DNA termini reached up to 117 pN during stretching. Notably, base pairs adjacent to intercalators were protected from strand separation, accompanied by additional helical unwinding. MM-PBSA/GBSA analysis revealed that the driving force for intercalation is the stacking energy, and the binding affinity was highest for minor groove binding. Persistence length decreased with single molecule binding but recovered with two molecules due to their electrostatic repulsion. Mechanical properties of intercalated DNA showed position-dependence, demonstrating that multiple intercalation modes coexist in solution. The heterogeneous nature of intercalation explains why experimental measurements reflect ensemble averages rather than single binding configurations.

Multi-objective reinforcement learning based on predicted structure feedback has been introduced into protein inverse folding. However, existing methods typically rely on a single model to optimize multiple structural objectives via a scalarized reward, which can bias the optimization toward dominant objectives and limit the exploration of diverse solutions. Here, we propose a online Symmetric Self-play Preference Optimization (SSP) framework that decouples the optimization of multiple structural objectives by training separate preference models with distinct reward signals, while enabling interaction through a shared sampling pool. This design allows the models to explore diverse optimization trajectories without enforcing a single dominant direction. Extensive experiments on both natural and de novo binder backbone inverse folding tasks demonstrate that SSP consistently improves sequence design self-consistency compared to single-model and existing baselines. Further analysis shows that different structural objectives are only partially aligned and induce distinct optimization directions, as evidenced by metric correlation and white-box analyses. This supports the effectiveness of decoupling objectives to enable higher design quality in protein design.

Bone morphogenetic protein receptor type IA (BMPR1A) is a key mediator of chondrogenesis and a validated therapeutic target for cartilage repair, yet existing BMP mimetic peptides suffer from low potency and the full-length protein (rhBMP-2) carries significant safety risks. Generative AI tools for protein design can now produce de novo peptide binders, but none have been applied to cartilage regeneration targets. Here, we benchmarked four architecturally distinct AI tools--RFdiffusion, BindCraft, PepMLM, and RFpeptides--to design candidate BMPR1A-binding peptides. We generated 192 candidates alongside 98 negative controls (290 total) and evaluated all complexes using AlphaFold 3 structure prediction, dual physics-based energy scoring (PyRosetta and FoldX), and contact recapitulation against the crystallographic BMP-2:BMPR1A interface (PDB: 1REW). A four-metric composite ranking identified a 15-residue PepMLM design (pepmlm_L15_0026) as the top candidate, combining favorable binding energy (PyRosetta dGseparated = -45.9 REU; FoldX {Delta}G = -19.4 kcal/mol) with the highest contact recapitulation among top-ranked peptides (11/30 gold-standard interface residues). Designed candidates significantly outperformed controls on ipTM (p = 0.002) and FoldX {Delta}G (p < 0.001). BindCraft candidates achieved the highest structural confidence (ipTM up to 0.81) but exhibited moderate contact recapitulation (mean 0.224), consistent with the computational hypothesis that they may engage alternative BMPR1A binding surfaces rather than the native BMP-2 interface. Physicochemical filtering yielded a shortlist of 54 candidates across all four tools. These results establish a reproducible computational framework for AI-guided peptide design targeting cartilage regeneration and identify specific candidates for future experimental validation via binding assays and chondrocyte differentiation studies. Author summaryDamaged cartilage has limited capacity to heal, and current biological therapies based on bone morphogenetic protein 2 (BMP-2) carry serious safety concerns including ectopic bone formation and inflammation. Short peptides that mimic BMP-2s interaction with its receptor BMPR1A could offer a safer, more targeted alternative, but designing such peptides from scratch is challenging. We used four different artificial intelligence tools--each employing a distinct computational strategy--to generate 192 candidate peptides designed to bind BMPR1A. We then evaluated all candidates using multiple independent computational methods to assess binding quality, energy favorability, and whether each peptide targets the correct site on the receptor. Our analysis identified a shortlist of 54 promising candidates, with a 15-residue peptide from the language model-based tool PepMLM emerging as the top-ranked design. We also found evidence that one tool (BindCraft) may produce peptides that bind BMPR1A at sites different from the natural BMP-2 interface, highlighting the importance of validating not just whether a peptide binds, but where it binds. Our computational framework and candidate peptides provide a foundation for future laboratory testing toward cartilage repair therapies.

Protein kinases are critical drug targets, requiring therapeutics that can modulate their active and inactive conformational states. While cofolding models can generate global folds directly from kinase sequences and ligand SMILES strings, these models have not yet been tested on their ability to recover ligand induced-fit conformational states of the kinase proteins. Here, we introduce KinConfBench, a curated benchmark of 2,225 high-quality human kinase chains to evaluate the ability of three state-of-the-art cofolding models--Boltz-2, Chai-1, and Protenix--to recover both canonical and rare conformational states. We show that geometric success metrics of a ligand pose in the active site does not correlate strongly with the correct kinase conformational state, motivating a new set of dynamical benchmarks for assessing cofolding models. While all three cofolding models achieve [~]65-75% prediction accuracy for kinase conformational classification, they exhibit severe mode collapse when performing multiple inferences, show negligible structural diversity in sampling induced-fit motions, and display a prevalent "apo-drift" in which all three cofolding models predominately predict the kinase to be in its ligand-free state. Our results highlight that capturing ligand-induced protein conformational diversity, not just geometric fit, is critical for next-generation structure-based drug discovery.

De novo protein design often produces thermostable proteins that denature above 100 {degrees}C, which complicates the analysis of their stability. Thermostable proteins can be unfolded by combined chemical and thermal denaturation followed by global analysis of multiple melting curves. Here, we have developed CheMelt, a new online tool for global analysis of unfolding data via an intuitive graphical user interface. We use nanoscale differential scanning fluorimetry followed by CheMelt data analysis to dissect the combined thermal and chemical denaturation of thirty-five de novo designed protein binders. Fifteen present sufficient fluorescence changes to extract thermodynamic parameters of unfolding. These de novo designed proteins have systematically lower {Delta}Cp and m-values than comparable natural proteins, which implies that they expose fewer hydrophobic residues upon unfolding. We show that a high thermostability of a designed protein does not necessarily imply a high equilibrium stability; and demonstrate the potential of CheMelt in dissecting thermodynamic properties for protein design and engineering.

Disease-related missense mutations reshape protein conformational energy landscapes, thereby altering biological function. However, mechanistically linking sequence variation to changes in conformational dynamics remains challenging for both experimental and computational approaches. Here, we introduce an internal-coordinate-based, essential-dynamics-refined elastic network model (ICed-ENM) that improves the physical fidelity of normal modes while capturing subtle mutation-induced side-chain effects and preserving computational efficiency. By constraining bond-length and bond-angle fluctuations and refining mode subspaces against experimentally observed collective motions, ICed-ENM provides a stable, structure-encoded description of intrinsic protein dynamics. Building on this framework, we developed a systematic mutation-scanning analysis that quantifies mutation impact as changes in vibrational entropy, providing a dynamic measure of mutation-induced redistribution within conformational energy landscapes. Validation against all-atom molecular dynamics simulations demonstrates that residues predicted as mutation hot spots induce substantial reshaping of free-energy landscapes, consistent with altered intrinsic conformational dynamics. Extending this analysis across a curated protein structure dataset reveals global patterns of mutation sensitivity across diverse structural and physicochemical contexts. Notably, these trends align with large-scale public mutation datasets, suggesting that our framework captures features relevant to pathogenic variation. Together, ICed-ENM and the associated mutation-scanning pipeline provide a scalable and mechanistically interpretable strategy to identify mutation-sensitive regions and substitutions, offering deeper insight into how sequence variation reshapes functional conformational landscapes.

Peptide-protein interactions are often transient and structurally elusive, necessitating computational approaches to identify both binding sites and peptide conformations. PatchMAN, one of the leading but computationally expensive biophysic-based global peptide-docking protocols, addresses this challenge by treating peptide docking as a protein-folding problem, using structural motifs from solved structures as templates that are subsequently refined using Rosetta FlexPepDock. Here we present PatchMAN2, which introduces 1) strategic fragment filtering and 2) local docking modes that focus sampling on relevant surfaces or known binding regions, thereby reducing the high computational cost of the original implementation due to extensive refinement of many non-productive low-quality fragments. Benchmarking shows that PatchMAN2 removes [~]30-70% of unnecessary fragments while preserving accuracy, substantially reducing runtime and improving the practical efficiency of peptide-protein docking.

1Long QT syndrome Type 2 (LQT2) is a genetic disorder caused by missense mutations in the KCNH2 gene that encodes the potassium channel KV11.1. Previous studies have shown that most KV11.1 missense mutations with loss-of-function phenotypes result from impaired trafficking from the endoplasmic reticulum to the plasma membrane. To investigate the molecular basis of these defects, we used molecular dynamics simulations to analyze two sets of disease-associated missense mutations: those that suppress and those that maintain normal channel trafficking. We focused initially on the conformational and dynamics differences between wild-type and several mutants of KV11.1 via molecular dynamics simulations when two K+ were placed in the selectivity filter (SF). Our study reveals that missense mutations in the S4 helix allosterically disrupt the selectivity filter, a critical determinant for proper channel trafficking. Trafficking-competent variants largely retained a wild-type selectivity filter structure, whereas trafficking-deficient mutants exhibited pronounced structural perturbations in this region. These findings suggest that certain LQT2-associated missense mutations in KCNH2 impair channel trafficking by compromising the structural integrity of the selectivity filter. We additionally found that second-site variants Y652C in the drug binding vestibule can correct structural defects associated with some mistrafficking variants.